source("intro_krig.r", echo=T)

library(rgdal)

library(gstat)

d <- readOGR("dtw_utm.shp", "dtw_utm")
OGR data source with driver: ESRI Shapefile 
Source: "dtw_utm.shp", layer: "dtw_utm"
with  612  rows and  12  columns

g <- readGDAL("srtm_utm_clip.tif")
srtm_utm_clip.tif has GDAL driver GTiff 
and has 80 rows and 81 columns

names(g) <- "elev.meters"

g$elev.srtm <- g$elev.meters * 3.2808399

d$elev.gw <- d$ELEV - d$DEPTH

d$elev.srtm <- overlay(g,d)$elev.srtm

d$elev.gw.srtm <- d$elev.srtm - d$DEPTH

d$elev.diff <- abs(d$elev.srtm - d$ELEV)

plot(d$elev.srtm, d$ELEV)

d$x <- coordinates(d)[,1]

d$y <- coordinates(d)[,2]

g$x <- coordinates(g)[,1]

g$y <- coordinates(g)[,2]

# krig an elevation surface ########################
v.elev <- variogram(ELEV ~ x+y, data=d)

plot(v.elev)

vm1.elev <- vgm(nugget=200, model="Exp", psill=3000, range=8000)

vm2.elev <- fit.variogram(v.elev, vm1.elev)

plot(v.elev, model=vm2.elev)

#uk.elev <- krige(formula=ELEV ~ x+y, locations=d, newdata=g, model=vm2.elev)

plt = spplot(uk.elev["var1.pred"], col.regions=bpy.colors(),
         scales=list(draw=TRUE, cex=0.7),
         sp.layout=list("sp.points", pch="+",

print(plt)

# krig a the groundwater  surface ########################
v <- variogram(elev.gw ~ x+y, data=d)

plot(v)

vm1 <- vgm(nugget=400, model="Exp", psill=2400, range=7000)

vm2 <- fit.variogram(v, vm1)

plot(v, model=vm2)

uk <- krige(formula=elev.gw ~ x+y, locations=d, newdata=g, model=vm2)
[using universal kriging]

plt = spplot(uk["var1.pred"], col.regions=bpy.colors(),
+          scales=list(draw=TRUE, cex=0.7),
+          sp.layout=list("sp.points", pch="+",  .... [TRUNCATED] 

print(plt)

library(rgl)

rgl.points(coordinates(d)[,1],d$elev.gw*10, coordinates(d)[,2], size=3,col="black")

cellsize =  summary(g)$grid$cellsize[1]

xdim = summary(g)$grid$cells.dim[1]

ydim = summary(g)$grid$cells.dim[2]

xorigin =  summary(g)$grid$cellcentre.offset[1]

yorigin =  summary(g)$grid$cellcentre.offset[2]

rgl.surface( ((1:xdim)*cellsize)+xorigin,
+              ((ydim:1)*cellsize)+yorigin,
+              uk$var1.pred*10, col="blue")

rgl.surface( ((1:xdim)*cellsize)+xorigin,
+              ((ydim:1)*cellsize)+yorigin,
+              uk.elev$var1.pred*10 + 3000, col="green")

